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ABSTRACT 

Robust modelling of strong lensing systems is fundamental to exploit the information they 
contain about the distribution of matter in galaxies and clusters. In this work, we present 
Lensed, a new code which performs forward parametric modelling of strong lenses. Lensed 
takes advantage of a massively parallel ray-tracing kernel to perform the necessary calcula¬ 
tions on a modem graphics processing unit (GPU). This makes the precise rendering of the 
background lensed sources much faster, and allows the simultaneous optimisation of tens of 
parameters for the selected model. With a single mn, the code is able to obtain the full poste¬ 
rior probability distribution for the lens light, the mass distribution and the background source 
at the same time. Lensed is first tested on mock images which reproduce realistic space-based 
observations of lensing systems. In this way, we show that it is able to recover unbiased es¬ 
timates of the lens parameters, even when the sources do not follow exactly the assumed 
model. Then, we apply it to a subsample of the SLAGS lenses, in order to demonstrate its use 
on real data. The results generally agree with the literature, and highlight the fiexibility and 
robustness of the algorithm. 

Key words: gravitational lensing: strong - methods: data analysis - methods: numerical - 
methods: statistical - techniques: image processing 


1 INTRODUCTION 

Strong gravitational lensing is now allowing researchers to address 
questions about the distribution of matter, both dark and light, 
within galaxy clusters and individual galaxies, that cannot be ad¬ 
dressed in any other way. The details of the distribution of mass in 
these systems have implications for such fundamental topics as the 
interactions of dark matter with itself and baryons, the mass of the 
dark matter particle and constraints on possible modifications of the 
theory of gravity, along with astrophysical questions such as how 
baryons segregate from dark matter and how the efficiency of star 
formation depends on the dark matter distribution. Strong lenses 
are also allowing astronomers to detect some of the most distant 
objects ever seen and study galaxy formation at the cosmic dawn. 
Future surveys such as Euclid (Laureijs et al. 2011) are expected to 
increase the number of known lenses by one or two orders of mag¬ 
nitude, and the scientific potential is indeed great if these lenses can 
be put to good use. 

At the heart of all strong lensing studies is the modelling 
of lenses based on astronomical images. For observations with 
resolved sources, Lefor, Futamase & Akhlaghi (2013) gave an 
overview of currently available software. Successful methods for 
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lens reconstruction have been presented recently by Warren & Dye 
(2003), Dye & Warren (2005), Suyu et al. (2006), Peng et al. 
(2006), Vegetti & Koopmans (2009), Newton et al. (2011), Ban- 
dara et al. (2013), and Nightingale & Dye (2015). 

Among the range of current lens reconstruction algorithms, 
forward reconstruction methods calculate an expected image from 
a model of the lens and source. By changing the parameters of the 
model, the forward reconstruction tries to minimise differences be¬ 
tween the reconstruction and observation, which are usually evalu¬ 
ated by a chi-square test or similar. These methods have been very 
successful in the reconstructions of recent galaxy-galaxy lensing 
surveys such as SLAGS (Bolton et al. 2006, 2008; Auger et al. 
2009). Because both the deflection field and the surface brightness 
distribution are modelled, these methods use all available informa¬ 
tion in the observation. A further advantage is the ability to use 
parametric models for background sources, thereby allowing a sys¬ 
tematic study of the parameters of lensed galaxies that would oth¬ 
erwise be inaccessible (Marshall et al. 2007; Newton et al. 2011; 
Bandara et al. 2013). The difficulty of the forward method lies in 
the thorough exploration of the model space, while simultaneously 
allowing for enough freedom in the model to represent reality. 

In this article, we present Lensed,^ a new implementation of 
the forward reconstruction method with focus on speed, quality of 
the lensing simulation, and statistical robustness of the results. The 
code was developed with the following features and goals in mind: 
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• User friendliness - We want the code to be widely and easily 
used by the scientific community. This includes easy importation 
of images, PSFs and masks. 

• Full posterior reconstruction - We want the code to report sta¬ 
tistically well-justified errors, including degeneracies, for the lens 
model parameters. 

• Insensitivity to initial parameters - We want to make it unnec¬ 
essary to guess the values of any parameters before beginning the 
fit. 

• Flexible modification of models - We want users to be able to 
easily create any model they choose and use the code to constrain 
its parameters. 

• Simultaneous fitting of sources and lenses - The number of 
sources and lenses can be chosen by the user. Sources include vis¬ 
ible foreground components and the sky brightness, both of which 
are usually subtracted in a separate step before lens fitting. 

• GPU acceleration - To achieve the necessary exploration of 
the parameter space in a reasonable period of time, the code must 
be fast. 

• Portability - To allow for wide use, the code must be able to 
run on different platforms with minimal reconfiguration. 

How we addressed each of these goals will be detailed below. 

The outline of this article is as follows. We describe the theory 
of the forward reconstruction method in section 2. The particulars 
of our implementation of the method are given in section 3. Sec¬ 
tion 4 contains the tests we perform on mock observations to verify 
and validate the performance of our code. Finally, in section 5 we 
use Lensed to analyse real space-based observations and compare 
the results to existing reconstructions. 


2 THE BAYESIAN FORWARD METHOD 


of light due to the lens, the surface brightness distribution Ig in the 
source plane is not directly observable. Instead, the surface bright¬ 
ness distribution 1 in the image plane is 

/(x) = 4(x) + /5(y(x)), (2) 


where y(x) must be evaluated according to the lens equation (1), 
and II is a contribution of (unlensed) light from the lens plane itself. 

Using the total surface brightness I(x) of the image plane, it is 
possible to calculate the expected luminosity Li, i - 1,..., m, for 
each of the m pixels of the image, as the integral 

L, = J I(x)dx (3) 

of the surface brightness I{x) over the pixel area P/. 

In any real observation, the obtained image is degraded by the 
instrument and the measurement process. This is usually modelled 
using a point-spread function (PSF), which is the observed image of 
an idealised point source. The PSF is usually provided for a given 
instrument and observation. In order to apply the PSF P, it is con¬ 
volved with the physical surface brightness distribution /, 


I\x) 




r) I{x + r) dr , 


(4) 


to give the effective surface brightness distribution F that is seen 
by the instrument. The expectation values for the pixel luminosity 
in this case become 


L^ = 


I r(x)dx , 

Jpij 


(5) 


which is of the same form as expression (3) without PSF. In the next 
section, we will not treat the reconstruction with and without PSF 
separately; it is understood that the convolution must be performed 
and the latter expression used whenever a PSF is to be applied. 


In this section, we wish to briefly develop the two fundamental 
parts of the forward method for lens reconstruction: simulation of 
the physical system - lenses, sources, foregrounds etc. - using an 
assumed model, and comparison of the predicted image with the 
true observation. The latter task requires careful attention, because 
for any comparison, one must first define a measure of similarity. 
Instead of using an ad hoc criterion, we seek an objective metric 
for our reconstructions, which is derived in section 2.2 in the form 
of a likelihood function for the model and data. With this, we can 
perform a reconstruction using the forward method in a fully prob¬ 
abilistic, Bayesian approach. 

2.1 Simulating the model 

The first part in modelling the physical system is to divide it into 
distinct planes on which there are two-dimensional lenses and 
sources. Lenses act as deflectors of the line of sight, and sources 
contribute to the plane’s surface brightness distribution. This is the 
usual approach to lensing, and an introduction can be found e.g. in 
the book by Schneider, Kochanek & Wambsganss (2006). 

We use a standard gravitational lens system consisting of one 
source plane and one image plane. The lens acts on the image plane; 
it maps an apparent angular position x to the angular position y on 
the source plane according to the lens equation 

yCx) = X - a{x) . (1) 

The function a is the deflection angle, which for the purpose of this 
method characterises the lens completely. Because of the deflection 


2.2 Comparing model and observation 

Having found the luminosity Li from (3) or (5) for each pixel, the 
expected values of the model can be compared to the array of ob¬ 
served pixel values d - {d\,..d^]- Because of the nature of the 
observational process, each individual di is a random variate. If the 
assumed model describes the physical process correctly, the mean 
of pixel value di is given by the expected luminosity Li in suitably 
chosen units, i.e. 

{di) = Li . (6) 

In addition to the values di, we require the variances 

= VarCd,) (7) 

for each pixel. These are sometimes provided with the observation, 
but often they are not and we must estimate the variance from the 
data itself as detailed in section 3.3. Then we can approximate the 
distribution of the pixel values by the normal distribution 

( 8 ) 

given that our model is indeed the correct one. Here we have made 
the implicit assumption that the ensemble of pixels is independently 
observed, meaning that the number of counts registered in one pixel 
does not depend on the counts for any other pixel. This assump¬ 
tion does not preclude the application of a PSF, which is an effect 
that happens before the counts for a given pixel are registered and 
is taken into account by using (5). However, correlations between 
pixels can be (and most certainly are) introduced by the instruments 
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or data reduction techniques such as drizzling (Fruchter & Hook 
2002). Our assumption therefore requires that these correlations are 
small, and that the variance for each pixel is enough to describe the 
statistics of the observation adequately. 

With the variance fixed, either externally or from the data, 
we are now able to give the probability distribution function for 
observing the pixel value dt given an expected luminosity L/, which 
is the Gaussian 


1 / Hdj-Ljf 


For the whole set of observed values d - {di,..., dm] and expected 
luminosities L = {Li,..., L^}, we then find the likelihood 


P{d I i) = n P(di I L,) 

i=l 


( 10 ) 


of observing data d given L by the model. As discussed before, the 
step from (9) to (10) requires the assumption that the pixels are 
independent. 

Finally, we are now able to build the parametric model which 
we want to reconstruct using likelihood (9). Such a model can be 
comprised of multiple individual lens and source objects, which as 
a whole describe the deflection and surface brightness in terms of 
the n parameters ^ = {^i,... For given values of we can 
use expression (3) or (5) to calculate the numerical values of the 
expected pixel luminosities Li in our image, which thus become 
implicitly functions of the parameters ^ of the model. Therefore, 
we can identify the likelihood of the parameter values and (10), 
and write 


m 


p{d I ^) = p(d IL) = ]q p{di I Ld. 

i=\ 


( 11 ) 


Using Bayes’ theorem, we subsequently invert the conditionality in 
likelihood (11) and find 


I d) 


P(d) 


( 12 ) 


for the posterior probability P(^ | d) of the parameters ^ given the 
observation d and prior probabilities P(^). The constant of normal¬ 
isation P(d) is the Bayesian evidence, which can be calculated as 
the marginal distribution 


p{d) = jp(d\^)p(^)d^ 


(13) 


and used for model selection. Because the evidence P{d) is a global 
constant, many algorithms that sample posterior (12) do not require 
its knowledge. This is in particular true for the class of Markov 
chain Monte Carlo (MCMC) samplers. 


3 IMPLEMENTATION 

We will now give an explanation of our implementation of the for¬ 
ward reconstruction method as detailed in section 2. A high-level 
algorithmic version could be the following: 

(i) Start with an image of observed pixel values dt and their ob¬ 
servational variance s]. If no variance is provided, generate it from 
the data. Apply mask when one is provided. 

(ii) Build a model of the lens and sources, parametrised by a 
number of parameters ^ with prior probabilities P(^). 

(iii) Pick a set of parameters ^ from the prior P(^). 


(iv) For the chosen parameters, calculate the expected luminosi¬ 
ties Li using (3). 

(v) Calculate likelihood P{d \ of the model using (9) and (10). 

(vi) Calculate the posterior probability P(^ | d) using the likeli¬ 
hood P{d I ^) and prior P(^). 

(vii) Repeat steps (iii) - (vi) until the parameter space is suffi¬ 
ciently sampled. 

Even though the algorithm in this form is conceptually simple, 
there has so far been no general purpose standard implementation 
of the forward reconstruction method for lensing. We think that 
there are two key obstacles that have to be overcome, which will be 
detailed in the following. 

The first problem is step (iv) of the schematic algorithm laid 
out above: A precise calculation of the expected luminosities (3) 
of the model requires a full ray-tracing simulation of the field of 
view to find the lensed surface brightness distribution, which must 
be subsequently integrated for each pixel of the image. With mod¬ 
ern space-based observations of galaxy-galaxy lensing events, the 
number of pixels that have to be calculated in this way routinely 
ranges from 10.000 to 100.000, making a numerical integration 
very costly. A possible way around this problem is the “poor man’s 
ray tracing” employed e.g. by Marshall et al. (2007), where the sim¬ 
ulated image samples the lensed surface brightness distribution at 
a fixed increased resolution. Naturally, there are limitations to the 
amount of time one can save in the computations before the method 
breaks down, as the real and estimated pixel luminosities near the 
highly magnified parts of the image depend critically on the quality 
of the numerical integration. 

The second problem of the forward method is step (vii), i.e. 
the sampling of the parameter space. For example, a reconstructed 
model might only fit some of the multiple images of a source at a 
time. In order to fit the remaining images, an algorithm must be able 
to move away from “good” parameter values, because the poste¬ 
rior distribution is naturally multi-modal. This rules out a standard 
maximum-likelihood approach for reconstruction, and requires a 
full sampling of the posterior instead. There are also a number of 
inherent degeneracies: for example, to swap the major and minor 
axis of a nearly-circular ellipsoidal lens, the position angle must be 
changed by +90°, even though the models are physically similar. 
Finally, the reconstruction is beset with strong correlations between 
parameters. Some are inherent in the surface brightness distribu¬ 
tion of sources, e.g. between the position, ellipticity, and position 
angle. Others arise from gravitational lensing, where observed fea¬ 
tures of an image can be attributed to both the light defiection by 
the lens and the intrinsic shape of the background source. Many 
strategies exist to improve the sampling of the parameter space, 
e.g. by taking known degeneracies into account. Further speed up 
could be achieved by a decomposition of the parameters into “fast” 
and “slow” variables, where “fast” variables are updated more fre¬ 
quently while the slow variables are kept fixed. This could help 
find e.g. the magnitude of sources quickly, as pixel values can be 
linearly scaled for any given values of the remaining parameters. 
Multi-modality can be tackled e.g. by first scanning the full pa¬ 
rameter space for global maximum-likelihood regions, which are 
subsequently analysed as usual (Birrer, Amara & Refregier 2015). 

With Lensed, we provide an implementation of the forward 
method that tries to handle the above issues in full generality. We 
believe that avoiding special treatment is the best way to achieve 
our stated goals of correctness of the simulation and robustness of 
the results. The result is an algorithm built around a massively par¬ 
allel ray-tracing kernel that performs the necessary calculations on 
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a modern graphics processing unit (GPU). Device code is written 
in OpenCL (a dialect of C), to ensure compatibility with commonly 
used devices and preserve the ability to compute on traditional CPU 
architectures when necessary or desirable (an overview of OpenCL 
was given by Stone, Gohara & Shi 2010). Our architecture has the 
added benefit that new source and lens models can be implemented 
quickly and without the need to recompile the code, opening up 
the way to purpose-built models for individual reconstructions. The 
setup of the physical system, as well as the program itself, can be 
done entirely with one configuration file. 


3.1 Source integration 

The simulation of a given model consists of the calculation of the 
expected pixel luminosities (3), which must be done numerically. 
Even though one cannot construct simple analytical models that 
capture the small-scale structure of real galaxies, it is nevertheless 
important that the integration is precise enough to not degrade the 
results of the reconstruction. Crude methods can introduce extra 
correlations between the source and lens components: For example, 
when using a single ray to determine the luminosity of a pixel, the 
lens might be coerced to optimally position the probe rays over 
the brightness distribution of the source. This can lead to a bias in 
results for an individual reconstruction, especially when a small or 
no PSF is used. We have found that the computational expense of 
performing reasonably precise integration of the pixel luminosities 
is often compensated by an increase in smoothness of the posterior 
distribution, which improves the overall sampling. 

The pixel luminosities (3) are thus integrated numerically by 
a weighted sum 

r ^ 

I{x)dx^YwkI{Xk). (14) 

with weights and abscissae x^ e Pi prescribed by some rule for 
numerical integration. In case a PSF is to be applied, it is assumed 
to be provided as an array of {2p + 1) x (2q -hi) function values 
{F_p ..., Fo,o,..., Fp^q] centred on the (0,0) element. The PSF 
is normalised automatically. After the pixel luminosities Lij (here 
indexed in two dimensions for simplicity) have been calculated, the 
integral (5) is replaced by the discrete convolution 

^0 = ( 15 ) 

k=-p l=-q 

This convolution is calculated directly on the device (GPU/CPU). 
Special care has to be taken near the borders of the image, where 
we pad the array to the required size by repeating the closest pixel 
value Lij. As usual, no significant signal should be present in the 
border region of p, q pixels when using a PSF. 

For the integration rule, we selected the Cartesian square of 
a classical seven-point GauB-Kronrod quadrature rule (Patterson 
1968) to be the default setting of Lensed. The choice is based on 
the testing procedure and comparison of different rules shown in 
section 4, and offers reasonable performance and good accuracy for 
the reconstruction of lenses as well as lensed background sources. 
As usual for GauB-Kronrod quadrature, an embedded three-point 
rule estimates the error of the integration for each pixel, which is 
saved with the output and can be used as a diagnostic for the quality 
of the numerical integration. 


3.2 Sampling of the parameter space 

For a practical application, it is almost always infeasible to explore 
posterior distribution (12) with a classic MCMC method, due to 
the high dimensionality of typical models and the usually strong 
correlations between individual model parameters. A more recent 
alternative to traditional MCMC samplers is the Nested Sampling 
algorithm given by Skilling (2004, 2006). We use the MultiNest 
library (Feroz & Hobson 2008; Feroz et al. 2009, 2013), which is 
an implementation and extension of the Nested Sampling algorithm 
that is well suited for working with the 10 < n < 50 parameters, 
multiple modes, and correlations that typically arise in strong tens¬ 
ing reconstructions. 

For a given problem, MultiNest requires the logarithm of the 
likelihood, which in our case (11) is 

logP(d I f) = Z 

Since the observational variance is fixed, we can use 
loglike(J I f) = -- V w,- (4- - Lif , (17) 

i=l 

to perform the required calculations, where we have introduced the 
weights Wi = 1 / 5 ^. We have dropped the second term from (16), 
which amounts to an overall normalisation that has no bearing on 
the results. An advantage of working with the weights Wf instead 
of the variance s] is that a pixel can easily be masked by setting 
its weight to zero. We further note that the sum in (17) amounts to 
the usual term, which we can use to summarise the ability of the 
given parameter values ^ to reconstruct the observed data. 

The MultiNest algorithm offers many settings that can be 
used to find the right balance between speed and thoroughness of 
the sampling for a given task, and Lensed exposes all of them to 
the user. We have tried to find reasonable default settings that work 
well in most typical situations, with a slightly larger emphasis on 
speed. These defaults should however in no way be considered to be 
ideal. We refer users to the relevant MultiNest documentation and 
invite them to find the combinations that work best for the specific 
problem at hand, as this often leads to significant improvements in 
speed, or quality, or - occasionally - both. 

We point out that in addition to the posterior distribution of 
the parameters, MultiNest calculates the evidence (13), which is 
a difficult problem in numerical Bayesian statistics. For any given 
model of lenses and sources, the evidence is a global, parameter- 
independent value that quantifies, as the name suggests, the confi¬ 
dence that the model is indeed responsible for the observed lensing 
event. It can thus be used to objectively compare different models 
that reconstruct the same observation. A good introduction into the 
topic was given by Kass & Raftery (1995). 

3.3 Variance estimation 

In the common case where the observation does not provide an 
independent estimate of the variance for each pixel, we must 
estimate it from the data df. When the lensing signal is sufficiently 
strong, we can estimate from the Poisson statistics of the photon 
flux. First, we compute the total photon counts ki for each pixel as 

ki = gi(di + bi) , (18) 

where gt is an effective gain and bt is an eventual offset that was 
subtracted from the input data. The effective gain is the conversion 
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factor from the data units to photon counts; for example, given an 
observation in ADU per unit time, the effective gain is 

gi - (electronic gain) X (exposure time) . (19) 

The offset bt is any constant that has been subtracted from the pixel 
values; this is often the estimated sky background. Assuming that 
the photons arrive at the detector with some fixed rate given by 
the physical system being observed, the counts kt follow a Poisson 
distribution, and we can estimate their variance as 

W2ix{ki)^ki^gi{di + bd . ( 20 ) 

Rewriting (18) for di and taking the variance, we find 

s] = VarW) = gf Var(^,.) = g7‘ (d,. + bd (21) 

as an estimator of the variance in case it was not independently 
provided by the observation. 

Where the signal from astronomical sources is low, other 
sources of noise, such as those due to the read-out process, become 
significant. They are not considered in the statistical variance de¬ 
scribed above and could be more difficult to estimate. Determining 
the noise level directly from the data, e.g. from the RMS values 
of empty patches of sky, depends highly on the data processing 
pipeline and can easily bias the reconstruction. Therefore, we do 
not include this option in the code, but advocate for careful ex¬ 
ternal preparation of the weight maps when necessary, so that the 
particulars of the observation can be correctly taken into account. 
Lensed still offers the possibility to include a multiplicative extra 
weight factor, for a quick modelling of other contributions to the 
noise. 


4 TESTING 

We have created a set of tests using virtual observations to verify 
and validate the Lensed code in a controlled environment. Each test 
contains a suite of 100 mock images of a galaxy-galaxy lensing 
event for a given lens and source model. While the parameters of 
the lens are fixed for all images, the source parameters are randomly 
selected for each image. By subsequently analysing all 100 mock 
images, we can separate the influence of the randomised source 
from the ability of Lensed to successfully reconstruct the fixed lens. 

In our tests, we consider a parametric lens and two different 
kinds of source model. In section 4.1, mock images containing 
sources with parametric profiles are used. Since both the lens and 
the source can be modelled exactly in Lensed, this provides an “up¬ 
per bound” on our reconstruction abilities, and we wish to verify 
the functionality of the code by recovering lens parameters with¬ 
out any significant bias. In the following section 4.2, we replace 
the parametric sources with more realistic ones based on observed 
galaxies. Because these sources cannot be modelled exactly, we use 
this test to validate the capability of Lensed to reconstruct the lens 
even if the sources are not described perfectly by their assumed 
model. The same data sets are used in sections 4.3 and 4.4 to in¬ 
vestigate the ability of Lensed to recover the parameters of sources. 
Finally, in section 4.5 we show how the choice of integration rule 
influences the results of these tests. 

In each test image, the lens is a singular isothermal ellip¬ 
soid (Kormann, Schneider & Bartelmann 1994) located at the 
centre of the image, with position xl - Jl - 0", Einstein ra¬ 
dius - 1.3624 arcsec, axis ratio qi - 0.75, and position an¬ 
gle 6l - 45° measured counter-clockwise between the major axis 
and x-axis. Assuming typical lens and source redshifts of zl - 0.2 



• 



# 

# 

# 
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Figure 1. Sample images from the mock catalogs of 100 strong lenses using 
parametric Sersic profiles (top) and observed galaxies decomposed into a 
shapelet basis (bottom) for the lensed galaxy. The lens in each is a singular 
isothermal ellipsoid centred on the 6 arcsec image, with Einstein radius 
ri = 1.3624 arcsec, axis ratio qi = 0.75, and position angle 6l = 45°. 
Background sources are positioned randomly within a disc circumscribing 
the caustic of the lens, with the remaining parameters randomised uniformly 
over their respective ranges. 

and = 1.0, respectively, the chosen Einstein radius translates 
into a velocity dispersion of cr^ = 250 km s"^ for the isothermal 
ellipsoid.^ 

A foreground galaxy is hosted in the lensing potential. Based 
on the majority of current galaxy-galaxy lensing observations, we 
assume it to be an early-type object, which is well parameterised 
by the profile of de Vaucouleurs (1948). The effective radius of the 
light profile is assumed to be r// = 2 arcsec, or about 1.5 Einstein 
radii, which is a reasonable ratio for typical lensing galaxies (e.g. 
Sonnenfeld et al. 2013). The other parameter values are randomly 
assigned for each image, using an axis ratio qu between 0.6 and 0.9, 
and major axis angle Oh between 25° and 65°, in order to maintain 
a certain alignment between light and mass profile. The foreground 
galaxy has a fixed magnitude of 16.5 mag and is co-centred with 
the mass distribution. 

The model of the background galaxy varies with the test we 
perform. In all cases, the background source is randomly positioned 
inside a circle circumscribing the caustic of the lens. This is done 
in order to observe the strong lensing features necessary for lens 
reconstruction. 

All mock observations in this section were realised with 
Glamer, a gravitational lensing simulator described by Metcalf & 
Petkova (2014) and Petkova, Metcalf & Giocoli (2014). Images are 
meant to realistically reproduce a space-based observation of the 
lensing systems. We use a HST-\ik& configuration for the telescope 
properties, simulating a 2000 second exposure in the F814W band. 
A typical PSF made with TinyTim^ is used for image synthesis 
and reconstruction. These settings match the observations analysed 
in section 5. The resulting images have a size of 120 x 120 pix¬ 
els at a resolution of 0.05 arcsec/px, for a side length of 6 arcsec. 
We choose to model this test on space-based observations because 
they are in general more difficult to reconstruct: Ground-based ob¬ 
servations of similarly-sized lenses contain fewer pixels, a larger 
PSF, and higher noise level, leading to a smoother likelihood with 
broader maxima. 

Fig. 1 shows a sample of the resulting virtual observations. 
Each test set of 100 mock images contains the commonly found 

^ However, Lensed works entirely in apparent units, and no information 
about cosmological distances is needed or implicitly employed. 

^ http://tinytim.stsci.edu/ 
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features of strong tensing such as rings, arcs, arclets, counter¬ 
images and crosses. While this is not supposed to be a statistically 
representative sample of the lenses observed in any specific survey, 
it should nevertheless cover a wide range of the configurations for 
which this code was intended, observed under realistic conditions. 


4.1 Lens parameters from an exact model 

For this test, the background galaxy follows an elliptical version 
of the Sersic profile (Sersic 1963; Graham & Driver 2005), 
defined by 

/(i?) = /o exp , ( 22 ) 

where R is the geometric mean radius 

R - ^ Jx ^ q + y ^ Iq (23) 

of the ellipse passing through the point x , y . The effective radius vs 
therefore describes the geometric mean of the semi-major and 
semi-minor axes of the half-light ellipse. The constant /q is cho¬ 
sen so that the source has a total magnitude of mag^. This profile 
is then translated to position V 5 , and rotated by the position an¬ 
gle 6s . Since our aim is not to create a realistic model, but rather one 
that can be recovered perfectly, the mock images contain only a sin¬ 
gle background component. Each image uses a different realisation 
of the uniformly random source parameters, with Sersic index ns 
between 0.5 and 8.0, effective radius vs between 0.1 and 0.4 arcsec, 
axis ratio qs between 0.1 and unity, and position angle 6s between 
0° and 180°. The range of values for the effective radius is com¬ 
patible with observational results reported by Shibuya, Ouchi & 
Harikane (2015), while the ranges for the other parameters allow 
for a very general representation of possible source configurations. 
The magnitude of the background galaxies before lensing is 23 mag 
in each case. 

Each mock image is separately analysed with Lensed using 
the default settings. In the spirit of this test, the model we assume 
is the same one used to generate the mock images: the lens is a 
singular isothermal ellipsoid, the foreground galaxy follows a de 
Vaucouleurs profile, and the background galaxy follows a Sersic 
profile. In total, the model contains 5 + 6 + 1 - 18 parameters. 
None of the parameters, including lens and source positions, are 
fixed or otherwise constrained. We use uniform prior distributions 
that span the whole range of possible parameter values. 

On a current workstation computer, the reconstruction of the 
full set of 100 images completes in a few hours. In all cases but one, 
the recovered maximum-likelihood lens parameters are very close 
to the input values. There was a sole outlier (Eig. 1, top row, first 
column), which returns a minimal yf^/dof value of 6.3, indicating 
that the sampling of the parameter space terminated before the true 
maximum-likelihood region was found. Eurther inspection of this 
system shows no particular oddities with respect to the total sample, 
and we therefore conclude that the premature termination is a nu¬ 
merical fluke. Since the reconstruction involves a certain amount of 
randomness, such behaviour is occasionally expected, and can usu¬ 
ally be detected by the high value ofy^^/dof. A subsequent analysis 
of this case with settings that result in a more thorough sampling 
of the parameter space (increased number of live points or reduced 
tolerance) aligns the results with the rest of the set. 

The distribution of the reconstructed lens parameters over the 
set of mock observations are shown in Eig. 2. No systematic bias 
in the outcome is apparent, and the results are tightly distributed 
around the true values. This is substantiated in Table 1, which con¬ 
tains the sample mean and standard deviation of the recovered lens 
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Figure 2. Posterior distributions of the lens parameters using an exact {red) 
and inexact (blue) model for the background sources, see text. Shown are 
the position XL,yL, Einstein radius axis ratio qi, and position angle 6l 
of a singular isothermal ellipsoid. Input parameter values are indicated by 
a vertical line (black). Results are marginalised over 100 HST-likQ mock 
observations with randomised sources. 


Table 1. Results of the reconstruction of mock images using an exact and 
inexact model for the background sources. Shown are the parameters for 
position XL, yL, Einstein radius axis ratio qi, and position angle 6l. The 
quoted values are the sample mean and standard deviation of the marginal 
distribution over 100 randomised realisations of the lens system, as shown 
in Fig. 2 


parameter 

input value 

exact source model 

inexact source model 

XL [arcsec] 

0 

-0.0000 + 0.0020 

-0.001+0.011 

yL [arcsec] 

0 

0.0005 + 0.0021 

-0.001 + 0.013 

rL [arcsec] 

1.36235 

1.36205 + 0.00080 

1.3619 + 0.0059 

qi 

0.75 

0.7499 ± 0.0026 

0.750 + 0.016 

Ol^ 

45 

44.97 ± 0.20 

45.2+1.7 


parameters over the set of all 100 mock images. The accuracy of 
the results, i.e. the vicinity of true value and average reconstructed 
value, is remarkably high. The precision of the results, i.e. the scat¬ 
ter of the recovered values about the mean, is good as well: posi¬ 
tion, Einstein radius, axis ratio, and orientation are constrained far 
below the pixel, per cent, and degree level, respectively. 

This basic test confirms that Lensed is able to recover the lens 
in a controlled setting without introducing systematic biases. Us¬ 
ing the true model for reconstruction, we are able to recover the 
lens parameters almost identically. In the next test, it will be shown 
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how the constraints on lens parameters get widened once the back¬ 
ground source model is no longer perfectly matched. 

4.2 Lens parameters from an inexact model 

A second set of images is made, in which background sources do 
not follow an analytical profile. Instead, they are extracted from 
a catalogue of observed galaxies in the Hubble Ultra Deep Field 
(Coe et al. 2006) and subsequently decomposed into the shapelet 
basis functions (Refregier 2003), using a code described by Mel¬ 
chior, Meneghetti & Bartelmann (2007). Shapelet functions form a 
complete and orthonormal basis which is suited for the extraction 
and the manipulation of galaxy images, and are described in detail 
by e.g. Meneghetti et al. (2008). 

Only galaxies with a redshift between 0.8 and 1.5 are selected 
from the Ultra Deep Field catalogue, to ensure a sample of galaxy 
shapes consistent with an assumed source redshift of = 1.0. 
Moreover, by making a cut in apparent magnitude at 27 mag, we 
filter out galaxies that were observed with low S/N. Each object is 
appropriately rescaled to have a (pre-lensing) magnitude of 23 mag 
and the same apparent size it would have had at redshift = 1.0. 

In this test, the lensed galaxies have a variety of shapes and 
substructures that cannot be modelled exactly. Our goal is to check 
in how far this mismatch infiuences our ability to recover the lens 
parameters. For the reconstruction, we use the same model as in 
section 4.1. As before, all priors are chosen to be uniform over the 
whole range of possible parameter values. 

The resulting parameter distributions of this test are shown in 
Fig. 2. The loss in precision due to the inexact nature of the recon¬ 
struction is clearly visible in the increased scatter of the recovered 
parameters. Still, the distribution of the results, while not exactly 
Gaussian, is remarkably symmetric and centred on the true values. 
As shown in Table 1, the accuracy for the whole sample of 100 re¬ 
constructions remains far below the respective scales of one pixel, 
per cent, or degree. Even though the standard deviation of the in¬ 
dividual parameters has increased five- to tenfold over the exact 
case, the precision is still excellent, with constraints on the posi¬ 
tion and Einstein radius at the intrapixel level, the axis ratio at the 
per cent scale, and the orientation correct to within two degrees. 
We have checked that the precision of the results can be increased 
using higher quality settings for the algorithm, although not to the 
same extent as in the exact case. 

With this second test, we have shown that Lensed retains the 
ability to recover our given lens parameters even when the back¬ 
ground galaxies cannot be modelled exactly. Shapelet sources are 
diffuse, usually with multiple peaks and troughs in the brightness 
distribution, and are generally only poorly fit by a single Sersic 
component. Nonetheless, we were able to use Lensed to constrain 
the parameters of the given SIE lens model to a very high degree. 
As we expect the reconstruction of real observations to be limited 
first and foremost by the model uncertainty due to an insufficient 
description of the mass distribution, we conclude that the effects 
of overly simplified sources are only secondary. We are thus con¬ 
fident that the reconstructions of real observations, such as those 
presented in section 5, are not limited by our implementation of the 
forward method. 


4.3 Source parameters from an exact model 

In many studies the gravitational lens is used as a serendipitous 
telescope through which a faint, distant high-redshift source can 



A0S [deg] 

Figure 3. Posterior distributions of the source parameters using an exact 
{red) and inexact {blue) model for the background sources. Shown are the 
errors of the reconstruction A = (recovered) - (input) for the recovered po¬ 
sition xs , ys , effective radius rg , magnitude mag^, Sersic index ng , axis 
ratio qg , and position angle 6g of the lensed sources, which are modelled 
in both cases with a Sersic profile. Input parameters for the inexact sources 
were found as explained in the text. Results were obtained using the Lensed 
algorithm on HST-Yike, mock observations. The results are marginalised 
over ensembles of observations with randomised sources, see text. 


Table 2. Constraining power of Lensed for the parameters of background 
sources. Shown are results for the errors A in the recovered values, which 
are expected to vanish. The corresponding distributions are shown in Fig. 3. 


parameter 

exact source model 

inexact source model 

Ax 5 [arcsec] 

-0.0000 + 0.0010 

-0.002 + 0.027 

Ay^ [arcsec] 

0.0001 +0.0011 

-0.002 + 0.029 

Ar^ [arcsec] 

0.003 + 0.019 

0.008 + 0.034 

Amag5 

0.003 + 0.033 

0.03 + 0.19 

Ang 

0.01 + 0.25 

0.07 + 0.25 

^qs 

0.0009 + 0.0085 

-0.005 + 0.046 

^0L [°] 

-1 + 14 

1 + 15 


be observed (Marshall et al. 2007; Newton et al. 2011; Bandara 
et al. 2013). For these cases, the reconstructions of the source be¬ 
comes as important as the reconstruction of the lens, and both must 
be performed simultaneously. While it was shown in the previous 
tests that the lens reconstruction is robust even when the sources 
are not perfectly represented by the source model, the same is not 
necessarily true for the reconstruction of source parameters under 
the influence of lensing, and the question merits further testing. 
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The first step is again validating that Lensed is able to recon¬ 
struct source parameters given an exact source model. For this, the 
same ensemble of mock observations are used. Here instead of us¬ 
ing the parameters ^ themselves, the error of the reconstructed pa¬ 
rameters, ^ - ^ 0 . is considered, so that all realisations have 
a common expectation value E[A^] = 0. The results are shown in 
Fig. 3, and a summary is given in Table 2. Photometric precision 
in this case of an exact model is 0.07 mag at the 2cr level, and the 
effective radius vs of the sources is constrained to 12% precision at 
the 2cr level. The least constrained parameter is the Sersic index ns , 
with a relative uncertainty of 10% at the 2cr level due to the gen¬ 
eral difficulty of recovering highly peaked sources (i.e. those with 
ns >2). Finally, the position angle 6s is very well reconstructed 
when the source is clearly elliptical, up to axis ratios of qs ~ 0.85, 
while for almost circular sources, the position angle is naturally 
unconstrained. 

Overall, Lensed can reliably reconstruct parametric sources 
which have been distorted by an unknown lens that is simultane¬ 
ously reconstructed. This is perhaps not a great surprise, but it is 
still a good validation of the fact that degeneracies and correlations 
between lens and source model do not restrict our ability of re¬ 
covering source parameters. The next step is to test whether this 
remains true when the sources are no longer described perfectly by 
the model used for their reconstruction. 


4.4 Source parameters from an inexact model 

When the true source is not perfectly described by the input model, 
we would still like to be confident that the reconstruction does not 
find systematically different source parameters due to some inter¬ 
action between the deflection and the source profile. However, in 
this case assessing the quality of the recovered source is more dif¬ 
ficult, because it is not obvious how to compute the input Sersic 
parameters of the sources. 

As described in section 4.2, the objects used in this test are real 
objects, which have been decomposed in shapelet basis functions. 
This means that they have an analytical profile which, in theory, 
extends out to infinity. It is not straightforward to define the best 
Sersic parameters for this kind of object. In order to perform the fit 
to a Sersic profile, it is necessary to produce a map of uncertain¬ 
ties that act as inverse weights in the fitting procedure. In practise, 
this is equivalent to producing mock observations of the objects, 
and analysing them with a reconstruction code such as the one de¬ 
scribed in this paper. We would like these Sersic ’’input parameters” 
to be similar to what one would get if the source were fit no lens, 
high signal-to-noise, no pixelisation and no PSF blurring. 

To this end, we produced virtual observations analogous to 
the ones presented in the beginning of this section for the lensing 
events. The only difference stands in the application of a uniform 
magnification that has the same total magnification |yu|tot as in the 
lensed case. The resulting images have high signal-to-noise val¬ 
ues, just as the lensed images we produced. Moreover, they are not 
pixelated or unresolved with respect to the point-spread function, 
because the original shapelet decompositions were derived from 
HST data, and thus they cannot contain any information on scales 
smaller than one of the original pixels. This scale, due the mag¬ 
nification applied, corresponds to ~ 10 pixels in the final images, 
significantly larger than the PSF smoothing kernel. Using the same 
total magnification as in the lensing case also gives us confidence 
that the same region of the galaxy will be typically dominant over 
the noise, and significant for the extraction of the Sersic parameters. 

Once the equivalent Sersic parameters of our sample of non- 
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Figure 4. Comparison of lens reconstruction results for different integration 
rules. Shown are the mean and standard deviations for the errors A over the 
sample of 100 reconstructions as a function of the chosen integration rule. 
For a list of the tested integration rules, as well as the parameters, see text. 


parametric shapelet sources have been found, the recovered and ex¬ 
pected parameters can be compared as in the case of parametric 
sources. The resulting marginal distributions are shown in Fig. 3. 
While the shape (position Xs, ys, axis ratio qs , position angle 6s ) 
and the luminosity (magnitude mag^) of the sources are less con¬ 
strained in comparison to the fully parametric case, the constraints 
on the parameters of the Sersic profile (effective radius rs and 
Sersic index ns) remain practically unchanged. This is substanti¬ 
ated in Table 2, which shows only a minor increase in variance for 
these parameters. 


4.5 Comparison of integration rules 

To conclude the set of tests, it is briefly investigated how the re¬ 
sults presented here are influenced by the choice of integration rule 
for the simulation of the model (cf. section 3.1). While it is impor¬ 
tant that errors in the calculation of the expected image are small 
enough not to degrade the results of the lens reconstruction, there is 
no need to waste time on a perfect integration when it can be shown 
that this does not improve the quality of the results. As usual in 
numerics, experimentation and testing is required to find the right 
choice for the task at hand, and an example is shown here. 

The first test of this section can be repeated using different in¬ 
tegration rules, in order to observe their influence on the recovered 
parameters. The results of such an experiment using the exact test 
set is shown in Fig. 4. The evaluated integration rules are 
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• point - a single point at the centre of each pixel, 

• sub2 - 2 X 2 subsampling of the pixels, 

• sub4 - 4 X 4 subsampling of the pixels, 

• g3k7 - GauB-Kronrod quadrature rule with 7x7 points, 

• gSkl 1 - GauB-Kronrod quadrature rule with 11x11 points, 

• g7kl5 - GauB-Kronrod quadrature rule with 15 X 15 points. 

As before, the model parameters for the lens are position 
scale radius ri, axis ratio qi and position angle PA/^, while the 
source has position xs,ys, effective radius rg, magnitude mag^, 
Sersic index ng , axis ratio qg , and position angle 6g. The results 
show a clear convergence of the results for the GauB-Kronrod rules, 
while 4x4 subsampling has similar accuracy when only the lens 
parameters are considered. Given the substantially lower number 
of integration points per pixel (16 instead of > 49), we conclude 
that 4x4 subsampling is the integration rule of choice for lens re¬ 
construction, at least for HST -like data using a realistic PSF at the 
same resolution. When source parameters are required, one of the 
GauB-Kronrod rules might be preferable in order to not bias the 
results, and we therefore select the 7 x 7-point rule as the default in 
our implementation. However, rule selection is flexible in Lensed, 
and we encourage users to experiment. 


5 APPLICATION TO REAL DATA 

In order to show the feasibility of the analysis of real images with 
Lensed, we consider a selection of lenses from the Sloan Lens 
ACS Survey (SLACS) catalogue (Bolton et al. 2006,2008). SLACS 
was a project that aimed to confirm with space-based observations 
lenses which were detected from ground-based SDSS spectroscopy. 
The final catalogue (Auger et al. 2010) comprises 85 confirmed 
lenses and 13 likely ones. Analysing the whole catalogue is out¬ 
side the scope of this work, and we limit ourselves to a demon¬ 
stration of the capabilities of Lensed on real observations. To this 
end, we analyse the first 8 observations of HST proposal 10886 in 
the F814W band that are available on the Hubble Legacy Archive. 
These SLACS lenses were originally investigated by Bolton et al. 
(2008), who modelled 7 of the systems, with the remaining system 
not modelled due to the presence of a companion to the lens galaxy 
(see Table 3). 

The base model for this analysis consists of a flat sky com¬ 
ponent (to account for imperfections in the process of sky subtrac¬ 
tion), a Sersic profile for the host galaxy, a SIE lens and a Sersic 
profile for the background galaxy. When the reconstruction with 
this model is not satisfactory (high we introduce additional 
Sersic components for the foreground or background galaxies. For 
this process, we consider a model detailed enough if there are no 
longer any systematic differences between the observation and re¬ 
construction, leading to just one or two components per visible ob¬ 
ject. In a more thorough analysis, one might choose the number of 
components using some objective criterion, and also impose a prior 
on their distribution (Brewer et al. 201 1). Since Bolton et al. (2008) 
forced the centre of the mass distribution to coincide with the cen¬ 
tre of the surface brightness distribution for all lenses, we restrict 
ourselves to the same assumption for the sake of comparison. We 
emphasise that this restriction is not necessary with our approach, 
as demonstrated in section 4. All other parameters are left uncon¬ 
strained. In those cases where we need to exclude additional visible 
objects close to the lens from the analysis, we provide Lensed with 
a suitable mask of the uninteresting regions. 

The estimated parameters of our reconstruction are given in 
Table 3, together with the number of source components used for 


foreground and background galaxies. Also listed are the results of 
the reconstructions by Bolton et al. (2008), which overall are in 
good agreement with our findings. The reconstructed images are 
shown in Fig. 5. As can be seen from the images, Lensed is able 
to model all individual lensing events, regardless of complexity, 
signal-to-noise of the background images and overlap between the 
different sources. A few cases required particular care due to their 
peculiar properties; they are briefly described in the following, in 
order to show the capabilities of the code and possible strategies for 
modelling non-trivial systems. 

J0808-I-4706 - This system was not modelled in the SLACS 
paper by Bolton et al. (2008) due to the presence of a close-by com¬ 
panion to the main lensing galaxy. It was qualitatively investigated 
by Orban de Xivry & Marshall (2009) as an example of a multi- 
component lens system, but without reconstruction. With Lensed, 
we used two separate and independent objects for the main lensing 
galaxy and its companion, modelling both as SIEs. We also used 
two Sersic components for the light of the main lensing galaxy, one 
Sersic component for the companion, and two Sersic components 
for the background source. In total, the model has 35 free param¬ 
eters, the highest number of all the tests presented in this work. 
Despite the complexity of the system, Lensed was able to explore 
the posterior (Eig. 6) and obtain reasonable results for the model. 
Of particular interest is the strong anti-correlation of the Einstein 
radii of the two lenses, a feature that would not have been evident 
in a maximum-likelihood reconstruction method. 

J0822-I-2652 - With the initial model, we were able to get a 
very accurate reconstruction of the bright arc on the right, but the 
residuals showed a clear sign of a circular structure around the cen¬ 
tral galaxy. As this could be the signature of a second background 
component, we added it to the model (as did Bolton et al. 2008). 
In the subsequent reconstruction, the evidence improved signifi¬ 
cantly (A log Ev ^ 800). Quoted are the results from this second 
model, although we cannot exclude the possibility that the residu¬ 
als in the first analysis were instead due to a mismatch between the 
light profile for the central galaxy and the model. If this is the case, 
the favoured model is much more elliptical, with an axis ratio of 
qL = 0.6285 + 0.0089. 

J0841-H3824 - In this case, most of the lensing information is 
contained in the very dim counter-image on the bottom left of the 
host galaxy. As explained below, we used a mask to make sure that 
Lensed makes use of it. We note that, even if in the reconstructed 
image looks similar to the one shown by Bolton et al. (2008), our 
recovered lens parameters differ significantly from theirs. This is 
probably due to the different method used to perform the delicate 
subtraction of the central galaxy. 

In the case of 10728-^3835 and 10841-^3824, the foreground 
galaxy is not well enough reconstructed by one or two Sersic com¬ 
ponents to reduce the significance of its residuals below that of the 
lensing signal. Therefore, with a completely free prior distribution 
for the position, the background source may be placed on top of the 
central galaxy to make up for the insufficiencies of the foreground 
model. This can be prevented by a better modelling of the fore¬ 
ground galaxy. However, a simple and efficient method is to use 
image plane priors for the background source positions. Instead of 
defining the prior distribution on the source plane, one can define a 
prior on the image plane for the position at which one of the images 
of the source is required to appear. Drawing from such a distribu¬ 
tion is done in a straightforward way by first sampling a position 
on the image plane, mapping it to the source plane, and assigning 
it to the source. This is a “prior” in the true sense of the definition, 
as it uses an observer’s prior knowledge regarding the apparent po- 
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Table 3. Estimated parameters for the SIE lens models of 8 SLACS lenses from HST proposal 10886. The table compares the mean and standard deviation 
obtained from Lensed with the SLACS lens modelling results published by Bolton et al. (2008). The columns Afg and Asrc give the number of components 
used to model the foreground and background galaxies, respectively. East-of-north angles from the original results have been changed to match our dehnition. 
The dotted fields indicate that the lens was not originally analysed. 


lens name 


Lensed 





SLACS 

q 



^siE [arcsec] 

q 

P.A. n 

Ng 

Asrc 

bsm [arcsec] 

P.A. r] 

Asrc 

SDSS J0029-0055 

0.9603 + 0.0010 

0.8824 + 0.0020 

113.09 + 0.45 

1 

1 

0.96 

0.89 

115.4 

2 

SDSS 10252-^0039 

1.0244 + 0.0006 

0.9226 + 0.0011 

16.07 + 0.47 

1 

1 

1.04 

0.93 

16.2 

3 

SDSS 10330-0020 

1.1012 + 0.0008 

0.8119 + 0.0032 

24.09 + 0.24 

1 

1 

1.10 

0.81 

23.2 

3 

SDSS 10728-^3835 

1.2344 + 0.0014 

0.8958 + 0.0027 

159.80 + 0.62 

2 

2 

1.25 

0.85 

157.6 

4 

SDSS 10808-F4706(a) 

1.1699 + 0.0022 

0.7991 + 0.0019 

9.37 + 0.41 

2 

1 





SDSS 10808-^4706(6) 

0.3335 + 0.0045 

0.9943 + 0.0024 

111.36+10.82 

1 

1 





SDSS 10822-^2652 

1.1081 +0.0037 

0.9577 + 0.0095 

99.14 + 7.39 

1 

2 

1.17 

0.88 

158.2 

2 

SDSS 10841-^3824 

1.3850 + 0.0024 

0.7831 +0.0019 

171.42+1.02 

2 

2 

1.41 

0.79 

1.4 

2 

SDSS 10903-^4116 

1.2911 +0.0022 

0.8784 + 0.0038 

68.04+ 1.08 

1 

1 

1.29 

0.90 

71.3 

2 
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Figure 5. Selection of SLACS observations from the F814W band and their reconstructions. Each reconstruction shows {left to right) the original observation, 
the maximum-likelihood model from Lensed, the lensed images of the background source only, and the residuals. Images are 5 arcsec by side with north up 
and east left. The scale ranges from -0.25 X {white) to X {black), where X is the 97th percentile image value of the recovered model for the first and second 
image, and the 99th percentile image value of the lensed background source for the third and fourth image (cf. Bolton et al. 2008). The overlays show masked 
regions {red) or image plane priors {green) if present in the reconstruction. 


sition of an image for the distribution of possible source positions. 
The use of image plane priors reduces the correlations between lens 
and source parameters, as a source normally has to be repositioned 
any time the lens is updated to make the observed images appear 
once again where they are observed. At the same time, the parame¬ 
ter space volume of the source position is dramatically reduced, as 


the observed position is well constrained from the start. The image 
plane prior regions we chose in the two cases are visible in Fig. 5. 

As evidenced by the reconstructions of this section, Lensed is 
able to obtain robust and meaningful results from the analysis of 
real data, where the uncertainties are much more severe than those 
of the mock data analysed in the previous section. The algorithm 
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Figure 6. Lens parameter constraints for the model of SDSS J0808+4706 (Fig. 5e). The two probable lenses are modelled as SIEs at the same redshift. The 
model contains a total of 35 free parameters, including the source and foreground surface brightness distributions. The triangle plot shows the marginalised 
one- and two-dimensional posterior distributions for the lens parameters, with contours at the 68%, 95%, and 99% confidence levels. 


is able to reconstruct at the same time the light coming from the 
central galaxy and from the background lensed sources. As in the 
analysis of mock data, no informative priors were used, and the 
analysis was, apart from the model choices described above, fully 
automatic. The few cases in which there are relevant discrepancies 
with previous analyses, these are due to inherent difficulties in the 
identification of lensing signatures, such as the circular residual in 
J0822-h2652 and the poor counter-image in 10841-^3824. These are 
peculiar cases that inevitably require some decision made by the 
modeller on the interpretation of the lensing event. In all cases in 
which there are clear signs of multiple images of the background 


source, Lensed is able to appropriately disentangle them from the 
light of the host galaxy, even in cases in which there is significant 
overlap (J0330-0020, J0903-f4116). It should be clarified that the 
formal errors quoted in this analysis are purely statistical, and cor¬ 
respond to real uncertainties only in the ideal case in which the 
model used for lenses and sources offers a correct description of 
these objects. In the analysis of real data, this may not be com¬ 
pletely true, and model-free quantities such as the Einstein radius 
of the lens can be affected by the uncertainty on the model. Thus, 
if one wants to interpret the quoted values in an extensive way, an 
additional empirical error of the order of a few percent should be 
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Table 4. Maximum-likelihood lens models for source reconstruction. The modelling procedure follows the one used by Bandara et al. (2013) as detailed in the 
text. The table shows maximum-likelihood parameters for the Einstein radius bsm, axis ratio q and position angle of the SIE model, as well as the magnitude 
7ext and position angle of the external shear. 


system 



Lensed 




Banda 

raetal. (2013) . 

q Text 


^siE [arcsec] 

q 

P.A. (SIE) 

Text 

P.A. (r) 

^src 

bsm [arcsec] 

-^src 

SDSS J0029-0055 

0.94 

0.84 

104.70° 

0.028 

168.90° 

2 

0.94 

0.86 

0.03 

2 

SDSS J0252-r0039 

1.02 

0.87 

27.71° 

0.029 

40.23° 

3 

1.03 

0.90 

0.02 

3 

SDSS J0330-0020 

1.10 

0.87 

11.02° 

0.022 

137.81° 

3 

1.10 

0.75 

0.03 

3 

SDSS J0728-r3835 

1.25 

0.72 

156.72° 

0.055 

152.91° 

3 

1.23 

0.74 

0.06 

3 

SDSS J0822-r2652 

1.14 

0.66 

155.90° 

0.087 

159.10° 

2 

1.12 

0.83 

0.05 

2 

SDSS J0841-r3824 

1.45 

0.66 

ijy 

0.010 

16.63° 

2 

1.40 

0.78 

0.01 

2 

SDSS J0903-r4116 

1.28 

0.91 

88.17° 

0.022 

137.47° 

3 

1.30 

0.87 

0.02 

3 


Table 5. Estimated parameters for the Sersic source models of 8 SLACS lenses from HST proposal 10886. The modelling procedure follows the one used 
by Bandara et al. (2013) as detailed in the text. The table shows the mean and standard deviation obtained from Lensed for the effective radius rhi, /-band 
AB magnitude m/ and Sersic index n of the background components. Repeated rows marked (a), (b), ... contain results for lens models with multiple source 
components. Bandara et al. (2013) quote cTmj ^ 0.38 mag and cr„ ^0.1 for the total systematic uncertainty of their magnitudes and Sersic indices, respectively. 





T .ENSED 


_ Bandara et al. f2013) 


system 

^src 

rhi [kpc] 

mi [AB] 

n 

rhi [kpc] 

mi [AB] 

n 

SDSS J0029-0055 

0.93 

2.688 + 0.152 

23.210 + 0.042 

2.728 + 0.100 

2.72 + 0.73 

23.51 

2.73 

SDSS J0252-r0039 

0.98 

0.814 + 0.014 

23.947 + 0.014 

1.487 + 0.039 

0.98 + 0.26 

24.43 

1.31 

SDSS J0330-0020 

1.07 

1.678 + 0.017 

21.981 +0.007 

2.585 + 0.035 

1.83 + 0.49 

22.00 

2.32 

SDSS J0728-r3835(a) 

0.69 

2.298 + 0.144 

23.680 + 0.057 

1.809 + 0.117 

2.87 + 0.77 

23.95 

1.99 

SDSS J0728-r3835(b) 

0.69 

7.923 + 0.565 

23.180 + 0.072 

3.963 + 0.123 

7.39+ 1.98 

24.31 

3.83 

SDSS J0728-r3835(c) 

0.69 

0.104 + 0.022 

27.671 + 0.104 

1.540+ 1.055 

0.11 + 0.03 

27.49 

0.04 

SDSS J0822-r2652(a) 

1.03 

1.163 + 0.032 

25.047 + 0.410 

0.361 + 0.041 

1.52 + 0.41 

25.36 

0.25 

SDSS J0822-r2652(b) 

0.59 

1.107 + 0.032 

27.789 + 0.037 

0.845 + 0.066 

2.01 + 0.54 

23.96 

0.72 

SDSS J0841-r3824 

0.59 

0.401 +0.013 

24.454 + 0.015 

0.061 + 0.007 

0.66 + 0.18 

24.43 

0.16 

SDSS J0903-r4116 

0.66 

1.716 + 0.015 

23.627 + 0.012 

0.547 + 0.017 

2.89 + 0.78 

23.51 

0.55 


considered, as done in other works (Bolton et al. 2008; Sonnenfeld 
et al. 2013). 

Finally, we wish to briefly comment on the reconstruction of 
background sources with Lensed. There is only limited room for 
comparison of our sample with the results of Newton et al. (2011), 
as only one of our systems has been analysed there. Instead, we 
look at the results of Bandara et al. (2013), where we And all our 
systems except J0808+4706 (for the same reason as above). To 
obtain comparable results, we perform a new analysis of our sys¬ 
tems following a similar reconstruction procedure. In a first step, 
we model the lens systems using a suitable number of foreground 
components, a SIE lens with external shear, and the same number of 
background components as Bandara et al. (2013). The maximum- 
likelihood parameters we obtain are shown in Table 4. We And that 
the results agree both with our previous analysis and the SLACS 
results quoted above. The differences in alignment and ellipticity 
are explained by the presence of external shear (see the discussion 
in Bandara et al. 2013). More importantly, our results are in good 
agreement with those obtained by Bandara et al. (2013), except for 
the position angles of the mass model and external shear. In light 
of the consistency of our results, we conclude that there is a prob¬ 
lem with their reported position angles,"^ which are therefore not 
reproduced here. Of note is the system J0822-I-2652, for which we 
And a more elliptical lens {q - 0.66) and a strong external shear 

^ If one (incorrectly) assumes that the individual images of our systems 
in the extended Eig. 3 of the online version of Bandara et al. (2013) are 
north-aligned, the position angles appear to be in agreement. 


(Text = 0.09). Considering that the signal of the second source was 
only identifled as a ring-like structure in the residual image (see 
above), it is likely that there is a pronounced degeneracy between 
ellipticity and external shear, and even subtle differences in the re¬ 
construction can lead to vastly different results in the degenerate 
q-y plane. This is corroborated by the close alignment of the SIE 
and external shear position angles. Besides, the level of agreement 
between the codes is generally high, particularly in light of the fact 
that our lens models are not built iteratively component after com¬ 
ponent as in Bandara et al. (2013), but instead in one pass. 

Once the lens systems have been modelled, Bandara et al. 
(2013) Ax the foreground and lens parameters to their maximum- 
likelihood values, and the individual background components are 
replaced by a single source. They perform this second step in order 
to obtain a “global” Sersic proflle model of the source plane. (This 
second step is not performed for systems with backgrounds clas- 
sifled as having “group” morphology.) Using Lensed, we are able 
to obtain these results by simply removing all but one background 
source components, except in the case of 10841-^3824. This sys¬ 
tem shows a strongly bimodal posterior distribution, in which the 
preferred solution reconstructs only one of the two sources. After 
suitably reducing the allowed parameter space, we obtain the de¬ 
sired “global” solution which envelopes both background source 
components. The full set of our results are shown in Table 5, where 
we once again And good agreement between the lens reconstruc¬ 
tion codes. Differences for the radii, magnitudes and Sersic indices 
of the reconstructed sources are reasonable, given the indicated er¬ 
rors. In the case of J0822-I-2652, we obtain smaller and dimmer 
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sources, which can be explained by differences in magnification 
for our more elliptical, strongly sheared lens model. Nevertheless 
we find acceptable agreement for the Sersic indices, and hence the 
shape of the source profile, in this system, as in the others. 

Based on this preliminary analysis, we conclude that Lensed 
is able to recover the properties of background sources. However, 
the size of our sample is too small to draw conclusions about 
whether there are intrinsic differences in the reconstruction of 
lensed sources when different algorithms are used. In the future, it 
will become necessary to perform a thorough and systematic com¬ 
parison of lensing codes for the analysis of lensed sources. 


6 CONCLUSION 

We have introduced Lensed, a new code for the forward lens re¬ 
construction from a parametric model. The forward method, pre¬ 
sented in section 2, offers conceptual simplicity and allows for a 
clear interpretation of the results for both lenses and sources. It is 
furthermore an important tool for testing new analytical models of 
lenses and - eventually - sources, at higher redshifts than otherwise 
accessible. 

There are however significant challenges that must be over¬ 
come in any truthful implementation. First, if the simulation of the 
assumed model is not accurate, i.e. an insufficient integration of the 
surface brightness distribution for each observed pixel, there will 
inevitably be a number of systematic biases in the results, as the 
quality of the simulation degrades close to the areas of high magni¬ 
fication, and thus high signal. Secondly, the parameter space for a 
lens reconstruction problem is usually highly degenerate with mul¬ 
tiple local maxima of possible solutions, making a straightforward 
sampling very difficult. 

In section 3, we have listed our solutions to the most severe of 
these problems. By carefully rewriting calculations in a massively 
parallel fashion, we were able to harness the multiprocessing capa¬ 
bilities of modern GPUs and maximise the efficiency of our code. 
The resulting increase in computational speed allowed us to simu¬ 
late the expected light distribution using numerical quadrature for 
each individual pixel, and still test thousands of individual param¬ 
eter settings per second. Coupling the fast and accurate simulation 
with a modern library for Bayesian analysis made it possible to 
achieve a true sampling of the system’s parameter space in reason¬ 
able time. 

Section 4 contains the tests we have performed to ensure our 
implementation is indeed working as expected. This has been done 
by analysing a fixed lens model in a sample of 100 mock images 
containing randomised sources. Recovering the lens parameters in 
a setting where the model can be exactly replicated in Lensed, 
we were able to show that no biases are introduced in the recon¬ 
struction. This remains true after switching to more complex back¬ 
ground galaxies based on observations, and thus an inexact model 
within Lensed. Even in this case, we recovered the true values to 
high accuracy. Thus we are confident that in real applications, the 
quality of the results is not limited by the algorithm we are using 
for the reconstruction. 

Finally, we used Lensed to analyse a number of real observa¬ 
tions in section 5. The data originates from the SLAGS survey of 
Bolton et al. (2006, 2008) and has been previously analysed there. 
We took this opportunity to compare our reconstructions to pub¬ 
lished results in the literature, highlighting eventual differences and 
their origins, while also showcasing how Lensed can be used in 
practise. 


There are a number of possible extensions to the algorithm 
in its current form. One such extension is support for the multi¬ 
plane lensing formalism of Petkova, Metcalf & Giocoli (2014), and 
Lensed was designed with this goal in mind. Another possibility is 
the simultaneous reconstruction of observations in multiple bands 
(Vika et al. 2013; Dye et al. 2014; Kostrzewa-Rutkowska et al. 
2014). Due to the speed and robustness of the implementation, it 
is also conceivable to use Lensed in a lens finding scenario such as 
the one recently proposed by Brault & Gavazzi (2015). Finally, it is 
possible to implement specialised treatment for lensing by galaxy 
clusters, the images of which are much larger in size than those 
of galaxy-galaxy lensing events, offering unique opportunities for 
optimisation. 
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